#delimit ;

/* SIMULATIONS */

/* We run the simulations for each type of distribution.  See the paper for each type of distribution.  */

quietly foreach type in 0 1 2 {;

  foreach b_true in 0 5 {;

    noisily display "`type'  `b_true'" ;

    drop _all;

/* We set the number of observations to ten thousand, which means that we repeat the simulations ten thousand times.  For each of these 10000 repetitions, we simulate 500 election outcomes */

    set obs 10000;

    if `type'==0 {;
      gen normalvote = rnormal(50,11) in 1/5000;
      replace normalvote = 15 if normalvote < 15 & normalvote != .;
      replace normalvote = 85 if normalvote > 85 & normalvote != .;
    };

    if `type'==1 {;
      gen normalvote = rbeta(2,5)*70 + 15 in 1/5000;
      replace normalvote = 15 if normalvote < 15 & normalvote != .;
      replace normalvote = 85 if normalvote > 85 & normalvote != .;
    };


    if `type'==2 {;
      gen     normalvote = rnormal(35,6) in 1/2500;
      replace normalvote = rnormal(65,6) in 2501/5000;
      replace normalvote = 15 if normalvote < 15 & normalvote != .;
      replace normalvote = 85 if normalvote > 85 & normalvote != .;
    };


    gen n = _n;

/* Here we define all of the underlying variables */
	
    foreach j in v  d  v_d1  v_abs  v_abs_d1  v_abs2  v_abs2_d1  v_abs3  v_abs3_d1  close  y  d_lo  close_lo  d_hi close_hi {;
      gen `j' = .;
    };

/* Here we define all of the underlying variables that are dependent on the bandwith used for estiamtions */

    foreach i in 1 2 3 4 5 40 {;
      foreach j in b_`i'  p_`i'  b_diff1_`i'  b_diff2_`i'  p_diff2_`i'  b_poly_`i'  p_poly_`i'  b_lin_`i'  p_lin_`i'  c_`i'  c_diff2_`i'  c_poly_`i'  c_lin_`i' {;
        gen `j' = .;
      };
    };
    gen b_optimal_bandwidth = .;
    gen c_optimal_bandwidth = .;
    gen   optimal_bandwidth = .;

/* Note, k defines each round of simulations */

    forvalues k = 1(1)10000 {;

      if mod(`k',1000)==0 {;
        noisily display `k';
      };

/* Here we reset all of the underlying variables */

      foreach j in v  d  v_d1  v_abs  v_abs_d1  v_abs2  v_abs2_d1  v_abs3  v_abs3_d1 {;
        replace `j' = .;
      };

/* Set the values from the simualtions to the underlying variables */

      replace v = normalvote + rnormal(0,9);
      replace d = (v > 50) if v != .;
      replace v_d1     = (v - 50) * (d==1);
      replace v_abs    = abs(v - 50);
      replace v_abs_d1 = abs(v - 50) * (d==1);
      replace v_abs2    = abs(v - 50)^2;
      replace v_abs2_d1 = abs(v - 50)^2 * (d==1);
      replace v_abs3    = abs(v - 50)^3;
      replace v_abs3_d1 = abs(v - 50)^3 * (d==1);

      replace y = normalvote + `b_true' * d + rnormal(0,9);

      foreach i in 1 2 3 4 5 40 {;

        local width = `i';

        foreach j in close close_lo d_lo close_hi d_hi {;
          replace `j' = .;
        };

        replace close = (v > 50 - `width' & v < 50 + `width');

        pause ;
  
/* Here we run all of the regressions and score the results for the line in the sample the correspond to k (i.e. the simulation number) */

/* BASELINE */

        reg y  d  if close==1;
        local b  =  _b[d];
        local se = _se[d];
        local df = e(df_r);
        local b_true_in_ci = ( `b' - `se' * invttail(`df', .025) < `b_true'  &  `b_true' < `b' + `se' * invttail(`df', .025) );
        test d = `b_true';
        local p = r(p);
        replace b_`i'      = `b'             if n==`k';
        replace p_`i'      = `p'             if n==`k';
        replace c_`i'      = `b_true_in_ci'  if n==`k';

/* POLYNOMIAL CONTROL FUNCTION */

        reg y  d  v_abs*  if close==1;
        local b  =  _b[d];
        local se = _se[d];
        local df = e(df_r);
        test d = `b_true';
        local p = r(p);
        local b_true_in_ci = ( `b' - `se' * invttail(`df', .025) < `b_true'  &  `b_true' < `b' + `se' * invttail(`df', .025) );
        replace b_poly_`i'  = `b'             if n==`k';
        replace p_poly_`i'  = `p'             if n==`k';
        replace c_poly_`i'  = `b_true_in_ci'  if n==`k';

/* LOCAL LINEAR CONTROL FUNCTION */

        reg y  d  v  v_d1  if close==1;
        local b  =  _b[d];
        local se = _se[d];
        local df = e(df_r);
        test d = `b_true';
        local p = r(p);
        local b_true_in_ci = ( `b' - `se' * invttail(`df', .025) < `b_true'  &  `b_true' < `b' + `se' * invttail(`df', .025) );
        replace b_lin_`i'  = `b'             if n==`k';
        replace p_lin_`i'  = `p'             if n==`k';
        replace c_lin_`i'  = `b_true_in_ci'  if n==`k';

/* PLACEBO CORRECTION */

        replace close    = (v > 50 - `width'/2 & v < 50 + `width'/2);
        replace close_lo = (v > 50 - `width' & v < 50);
        replace d_lo     = (v > 50 - `width'/2) if v != .;
        replace close_hi = (v > 50 & v < 50 + `width');
        replace d_hi     = (v > 50 + `width'/2) if v != .;

        reg y  d_lo  if close_lo==1;
        local b_lo = _b[d_lo];
        reg y  d_hi  if close_hi==1;
        local b_hi = _b[d_hi];
        local b_diff  = `b' - (`b_hi' + `b_lo')/2;
        replace b_diff1_`i' = `b_diff' if n==`k';

        reg y  d  d_hi  d_lo  if (close==1 | close_lo==1 | close_hi==1);

        local b = _b[d];
        local b_hi = _b[d_hi];
        local b_lo = _b[d_lo];
        local df = e(df_r);

        lincom d - (d_lo + d_hi)/2;
        local b_diff  = r(estimate);
        local se_diff = r(se);

        test d - (d_hi + d_lo)/2 = `b_true';
        local p_diff = r(p);
        local b_diff_x  = `b' - (`b_hi' + `b_lo')/2;
        local b_true_in_ci = ( `b_diff' - `se_diff' * invttail(`df', .025) < `b_true'  &  `b_true' < `b_diff' + `se_diff' * invttail(`df', .025) );
        replace b_diff2_`i' = `b_diff'       if n==`k';
        replace p_diff2_`i' = `p_diff'       if n==`k';
        replace c_diff2_`i' = `b_true_in_ci' if n==`k';

      };

      rdob y  v, c(50);
      local b  = $t_rd;
      local se = $rdse;
      local h_opt = $h_opt;
      count if v > 50 - `h_opt' & v < 50 + `h_opt';
      local df = _result(1) - 4;
      local b_true_in_ci = ( `b' - `se' * invttail(`df', .025) < `b_true'  &  `b_true' < `b' + `se' * invttail(`df', .025) );
      replace b_optimal_bandwidth = `b'            if n==`k';
      replace c_optimal_bandwidth = `b_true_in_ci' if n==`k';
      replace optimal_bandwidth   = `h_opt'        if n==`k';

    };


    compress;
    drop normalvote - close_hi;

/* We save the output from each distribution type individually */

    if `type'==0 {;
      save tmp_simulations_normal_`b_true', replace;
    };

    if `type'==1 {;
      save tmp_simulations_skewed_`b_true', replace;
    };

    if `type'==2 {;
      save tmp_simulations_bimodal_`b_true', replace;
    };

  };

};


